R Notebook

Loading Necessary Packages

Reading in Dataframes

# Community data wide format
cover.wide <- read.csv("../../Data/cover_wide_tradeoffs.csv")
# Site attributes
site.attr <- read.csv("../../Data/comb-by-plot-clim-soil-diversity-02-Aug-2019.csv")

Subsetting datasets to attribute and community matrices

Fitting RRPP models to each site

output_list <- list()
counter <- 1

sitenames = unique(attr.mat$site_code)
sitenames = sitenames[!(sitenames %in% c("doane.us", "pape.de", "sevi.us"))]

# For all unique sites selected
for(sites in sitenames){

  print(paste(counter,":", sites))

  # Subset to a single site
  com.subset <- cover.mat[attr.mat$site_code == sites,]

  # Remove all zero columns
  com.subset <- com.subset[,colSums(com.subset) > 0]

  # Pull out relevant plot attributes for each site
  attr.subset <- data.frame(attr.mat) %>%
    filter(site_code == sites) %>%
    mutate(year_trt = as.factor(year_trt),
           plot = as.factor(plot),
           block = as.factor(block))

  # Check that these two matrices are the same size
  expect_true(nrow(attr.subset) == nrow(com.subset))
  expect_true(nrow(attr.subset) > 0)

  mod_rrpp <- lm.rrpp(com.subset ~ year_trt + trt_K_num + trt_P_num + trt_N_num + block,
                      data = attr.subset,
                      SS.type = "III",
                      iter = num_iter,
                      print.progress = FALSE)

  # Saving ouput
  output = list(siteinfo = attr.subset,
                aovtable = anova(mod_rrpp, effect.type = "F",
                                  error = c("Residuals", "Residuals", "Residuals",
                                            "Residuals", "Residuals")),
                specscores = coef(mod_rrpp))

  output_list[[counter]] = output
  counter = counter + 1
}

storagelist <- list()
for(i in 1:length(output_list)){
  anovatab = bind_cols(site = as.character(unique(output_list[[i]]$siteinfo$site_code)),
            data.frame(matrix(output_list[[i]]$aovtable$table$`Pr(>F)`, nrow = 1)))
  colnames(anovatab)[2:ncol(anovatab)] =   rownames(output_list[[i]]$aovtable$table)
  storagelist[[i]] = anovatab
}
write.csv(x = bind_rows(storagelist),
  "../../Data/tradeoffs_RRPP_Pvals.csv")
storagelist <- list()
for(i in 1:length(output_list)){
  sp_scores = data.frame(output_list[[i]]$specscores[grep("trt_", rownames(output_list[[i]]$specscores)),])
  tomatch = c(4:10)
  # sp_scores = sp_scores[grepl(paste(tomatch, collapse = "|"), rownames(sp_scores)),]
  spectab = bind_cols(site = rep(unique(output_list[[i]]$siteinfo$site_code), nrow(sp_scores)),
                    trt = rownames(sp_scores),
                    sp_scores)
  storagelist[[i]] = spectab %>% mutate(trt = gsub("_num", "", trt))
}
specscores_full = data.frame(bind_rows(storagelist))
specscores_full[is.na(specscores_full)] = 0
write.csv(specscores_full, "../../Data/tradeoffs_specscores.csv")

Calculate dot-product of pairwise correlations

Given these estimated parameters, we can calculate the correlation between species responses to different pairs of nutrient treatments through the dot-product of the fitted vectors between treatments within a site.

After calculating this dot-product, the distributions of each set correlation scores can be seen in histograms below:

• Values close to zero indicate that species responses to two treatments within a site are not correlated with one another.

• Values close to one or negative one indicate that species responses to two treatments within a site are highly positively or negatively correlated.

# Normalization function -- squared sums of vector = 1
scalar1 <- function(x) {x / sqrt(sum(x^2))}

# Defining dot-product function
dotprod <- function(x, normalize = TRUE){

  # First, select the appropriate treatment and remove columns that aren't plant response
  N_vec <- x %>% filter(trt == "trt_N") %>% select(-trt, -site)
  P_vec <- x %>% filter(trt == "trt_P") %>% select(-trt, -site)
  K_vec <- x %>% filter(trt == "trt_K") %>% select(-trt, -site)

  N_change <- sqrt(sum(N_vec^2))
  P_change <- sqrt(sum(P_vec^2))
  K_change <- sqrt(sum(K_vec^2))

  # If normalizing, run the normalize function on each vector
  if(normalize == TRUE){
    N_vec <- scalar1(N_vec)
    P_vec <- scalar1(P_vec)
    K_vec <- scalar1(K_vec)
  }

  # Generate output dataframe of dot products (pracma::dot function)
  output = data.frame(NP = pracma::dot(as.numeric(N_vec), as.numeric(P_vec)),
                      NK = pracma::dot(as.numeric(N_vec), as.numeric(K_vec)),
                      PK = pracma::dot(as.numeric(P_vec), as.numeric(K_vec)),
                      N_change = N_change,
                      P_change = P_change,
                      K_change = K_change)

  return(output)

}

# Run over our dataset
dotoutput = specscores_full %>% group_by(site) %>%
  do(dotprod(.))

# Load permANOVA p-value output
site_pvals = read.csv("../../Data/tradeoffs_RRPP_Pvals.csv", header = TRUE, stringsAsFactors = FALSE)

# Join all together (sites, treatments, dot-product pairs, perMANOVA p-values)
dot_full = left_join(site_pvals %>% select(-Total, -Residuals), dotoutput) %>%
  select(-X) %>%
  rename("trt_K" = "trt_K_num", "trt_P" =  "trt_P_num", "trt_N" = "trt_N_num")

# Write final data product
write.csv(x = dot_full, "../../Data/dot_full.csv")

Plotting correlation between different vectors of change

dot_full <- read.csv("../../Data/dot_full.csv")
p1 = dot_full %>% # filter(trt_N < .05 & trt_P < .05) %>%
  ggplot(aes(x = NP,
             fill = trt_N < 0.05 & trt_P < .05)) +
  geom_histogram(aes(y = ..count..), color = "black", alpha = .5, position = "identity", bins = 25,
                 data = dot_full %>% filter(trt_N < .05 & trt_P < .05)) +
  geom_histogram(aes(y = -..count..), color = "black", alpha = .5, position = "identity", bins = 25,
                 data = dot_full %>% filter(!(trt_N < .05 & trt_P < .05))) +
  xlim(-1, 1) +
  ggtitle("N - P Correlation") +
  guides(fill = FALSE) +
  scale_y_continuous(limits = c(-4,4), breaks = seq(-6, 6, by = 2), labels = abs( seq(-6, 6, by = 2))) +
  geom_text(aes(x = -.25, y = 3.5), label = "Both vectors significant (p < 0.05)", color = "turquoise4") +
  geom_text(aes(x = -.25, y = -3.5), label = "=< 1 vectors significant (p < 0.05)", color = "tomato2")
p2 = dot_full %>% # filter(trt_N < .05 & trt_K < .05) %>%
  ggplot(aes(x = NK,
             fill = trt_N < 0.05 & trt_K < .05)) +
  geom_histogram(aes(y = ..count..), color = "black", alpha = .5, position = "identity", bins = 25,
                 data = dot_full %>% filter(trt_K < .05 & trt_N < .05)) +
  geom_histogram(aes(y = -..count..), color = "black", alpha = .5, position = "identity", bins = 25, 
                 data = dot_full %>% filter(!(trt_K < .05 & trt_N < .05))) +  
  xlim(-1, 1) +
  ggtitle("N - K Correlation")+
  guides(fill = FALSE) +
  scale_y_continuous(limits = c(-4,4), breaks = seq(-6, 6, by = 2), labels = abs( seq(-6, 6, by = 2)))
p3 = dot_full %>% # filter(trt_K < .05 & trt_P < .05) %>%
  ggplot(aes(x = PK,
             fill = trt_K < 0.05 & trt_P < .05))+
  geom_histogram(aes(y = ..count..), color = "black", alpha = .5, position = "identity", bins = 25,
                 data = dot_full %>% filter(trt_K < .05 & trt_P < .05)) +
  geom_histogram(aes(y = -..count..), color = "black", alpha = .5, position = "identity", bins = 25,
                 data = dot_full %>% filter(!(trt_K < .05 & trt_P < .05))) +
  xlim(-1, 1) +
  ggtitle("P - K Correlation")+
  guides(fill = FALSE)+
  scale_y_continuous(limits = c(-4, 4), breaks = seq(-6, 6, by = 2), labels = abs( seq(-6, 6, by = 2)))
grid.arrange(p1,p2,p3,
             top = textGrob("Site Response Correlations",gp=gpar(fontsize=20,font=3)),
             nrow = 3)

Visualizing Example Communities

rdaFit = function(site){
  com.subset <- data.frame(cover.mat)[attr.mat$site_code == site,]
  
  # Remove all zero columns
  com.subset = com.subset[,colSums(com.subset) > 0]
  
  # Pull out relevant plot attributes for each site
  attr.subset <- data.frame(attr.mat) %>% filter(site_code == site)
  
  # Check that these two matrices are the same size
  expect_true(nrow(attr.subset) == nrow(com.subset))
  
  # Generate and return an NMDS figure
  mod_out <- vegan::metaMDS(com.subset, distance = "euclidean", autotransform = FALSE, try = 50)
  env <- envfit(mod_out, attr.subset %>% select(block, trt_N_num, trt_P_num, trt_K_num) %>% mutate(block = as.factor(block)))

  # Generating GGplot
  ggplot() + 
    geom_point(aes(x = MDS1, y = MDS2), pch = 21, fill = "white", alpha = .5, size = 2,
               data = data.frame(mod_out$points) * attr(mod_out$points, "internalscaling")) +
    geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), data = data.frame(env$vectors$arrows), size = 2, 
                 arrow = arrow(length = unit(0.5, "cm"), type = "closed")) +
     geom_point(aes(x = NMDS1, y = NMDS2), data = data.frame(env$vectors$arrows), alpha = 0, size = 6) + 
    ggrepel::geom_text_repel(aes(x = NMDS1, y = NMDS2), 
              data = data.frame(env$vectors$arrows), 
              label = gsub("_num", "", rownames(data.frame(env$vectors$arrows))),
              size = 6) +
    ggtitle(site)
}

rdaFit("rook.uk"); ggsave("../../Figures/rook_nmds.pdf", height = 5, width = 5)
rdaFit("bnch.us"); ggsave("../../Figures/bnch_nmds.pdf", height = 5, width = 5)

Summarising communities by functional group abundance